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We develop a graphical method for computing the virial expansion coefficients for a nonrelativistic 
quantum field theory. As an example we compute the third virial coefficient &3 for unitary fermions, 
a nonperturbative system. By calculating several graphs and performing an extrapolation, we arrive 
at 63 = —0.2930, within 0.7% of a recent computation 63 = —0.29095295 by Liu, Hu and Drummond 
[1], which involved summing 10,000 energy levels for three unitary fermions in a harmonic trap. 
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I. INTRODUCTION 

The virial expansion allows one to express the equation 
of state of a non-ideal gas in a density expansion, and is 
equivalent to a fugacity expansion of the grand potential 
density at nonzero chemical potential: 

-^^PP^^[z + hz^ + hz' + ...] (1) 

or dn/dfi = -2(y/A3)^„n6„z". Here V, P and ^ 
are the volume, pressure and inverse temperature respec- 
tively, z — e^^ is the fugacity, and A = y/2T:(3/M is 
the thermal wavelength; the 6„ are dimensionless quan- 
tities directly related to the virial coefficients. The 0{z) 
contribution is independent of interactions, and therefore 
the ideal gas term has been factored out front (we 
assume here a gas of spin 1/2 fermions). The thermal 
wavelength provides a natural length scale, and the fu- 
gacity expansion is expected to be valid when A is short 
compared to the average inter-particle distance, and long 
compared to the range of interactions. This expansion 
is of current interest because of experimental focus on 
the properties of dilute atomic gases; and the case of 
fermionic atoms at a Feshbach resonance — where the 
two-body scattering length diverges — is of particular 
interest to both theorists and experimentalists. Such 
examples of "unitary fermions" are strongly interacting 
conformal systems interpolating between the BCS and 
BEC regimes, with universal properties that serve (on a 
completely different length scale) as an interesting start- 
ing point for effective field theory treatments of interact- 
ing nucleons [U |3] . There has been extensive theoretical 
interest in computing the parameter 63 in the fugacity 
expansion for unitary fermions [4j [5], culminating in a 
high accuracy determination 63 = —0.29095295 from a 
spectral study of the 3-fermion system, solving for the 
lowest 10,000 energy levels for three unitary fermions in 
a harmonic trap [T], a value that appears to agree with 
experimental results [7] . 
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It would be convenient to have a method for computing 
the bn coefficients directly using field theory techniques, 
particularly if accurate results could be obtained by com- 
puting a small set of diagrams. Despite the fact that H. 
is directly related to the sum of one particle irreducible 
Feynnian diagrams in a finite temperature field theory, 
the theory is not ideally suited to this task for several 
reasons: interparticle potentials are not in general read- 
ily described by Feynman rules; the sum of multiple in- 
terparticle interactions cannot typically be computed an- 
alytically or numerically without resorting to solving the 
corresponding Schrodinger or Lippman-Schwinger equa- 
tion; the Feynman graphs are functions of arbitrary fi 
and there is no simplification gained by the expansion in 
z. 

In this Letter we devise a graphical expansion that 
circumvents these difficulties, and demonstrate its util- 
ity by performing analytic calculations of 63 for unitary 
fermions; it has some features in common with the ap- 
proach of refs. [HIS]. The calculation is mostly analyt- 
ical, with some integrals performed numerically, and we 
will show that an extremely accurate determination of 
63 can be obtained. The two components of our general 
procedure are (i) to perform the "dual" of the Matsubara 
sum over discrete frequencies — in the sense of a Poisson 
resummation — which directly leads to a fugacity expan- 
sion; (ii) to use a dimer field as in ^8 , which is designed 
to reproduce the continuum 2-body phaseshift, thereby 
bypassing discussion of potentials and leading to purely 
local interactions in space. We consider these two inno- 
vations in turn, first addressing the case of free fermions, 
then including two-body interactions. 

II. CHRONOGRAPHS 

We consider a dilute gas comprised of a single species of 
nonrelativistic spin half fermion; generalization to bosons 
or more species is straight forward, but we have not con- 
sidered the relativistic case. In the Euclidian time formu- 
lation of finite temperature field theory, is given by sum 
of IPI vacuum Feynman diagrams, where the theory is 
analyzed for Euclidian time r compactified with period 
/3 and antiperiodic (periodic) boundary conditions im- 
posed for fermions (bosons). For a free spin ^ fermion, 
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FIG. 1: The free fermion contribution to fi. The conventional 
finite temperature Feynman diagram on the left is expanded 
as a sum over worldline loops about the compact time di- 
rection ("chronographs") with winding number v making a 
contribution proportional to z" . The black dot indicates the 
nontrivial topology, and Euclidian time increases in the coun- 
terclockwise direction. 



FIG. 2: The contour C in eq. ([5| is designed to pick up 
contributions from all cuts and poles along the real energy 
axis. 



is given by the one loop diagram on the left in Fig. [T] It 
is convenient instead to compute d^/d^ with the result 

5^--^TrG„ --— y — ^GoK,p) (2) 

where G^(a;,p) = e*'^°^/(iaj — {e^ — ^)) is the free Eu- 
cHdian propagator, Cp = p^/2M, and aj„ = 27r(n-|- 1)//3; 
the factor of 2 is from the two spin states, and the mi- 
nus sign from the fermion loop. The frequency sum is 
trivial to compute and the result can be subsequently 
expanded in powers of the fugacity, but it is interesting 
to note that a Poisson resummation yields the fugacity 
expansion directly (see also appendix B of Ref. [S]): 

^^Go^(c.„,p) = ^(-irG^(i./3,p) 

n u 

oo 

= ^{-'^)^'"^^^z''e'"'^'^ (3) 

where G^(r, p) = —9{t — 0+)e^'^(^p^^^ is the Fourier 
transform of G^(w,p). The Poisson formula exchanges 
the sum over Matsubara frequencies for a sum over the 
winding number v of worldlines wrapping around the 
compact time direction, each term proportional to z'^ , as 
shown graphically in Fig. [l] We therefore immediately 
read off the 6„ coefficients for a free fermion: 

^(1) _ i_i\n+^^ f e-"'^^p - '^''^^"^^ (4) 

" ^ ' n J (27r)3 " n5/2 " 

Note this result includes both the (—1) from the Feynman 
graph, as well as a factor of (— 1)"^ from fermion worldline 
loops due to antiperiodic boundary conditions. 

We will refer to the diagrams on the right in Fig. [T] 
as "chronographs" , which allow one to compute directly 
the n*'' term in the fugacity expansion of fl or dVl/dfi. 
The rules for chronographs can be easily generalized for 
computing the fugacity expansion in interacting systems: 
(i) Chronograph propagators Q (r, p) can be defined in 
terms for the Minkowski propagator G^(i?,p) =i/{E — 
Ep -\- ie) — Q via the contour integral along the path 
C shown in Fig|2j which simply picks up all the physical 
poles and cuts: 

g(r,p)^0(r-O+) |^ge-^^G^^(i?,p) (5) 



Q should be thought of as a multi-valued function of t 
living on a compact manifold of circumference /3. For a 
free fermion, Qi^{t,'p) = — 6'(T)e~'^'^p; (ii) vertices (from 
the Euclidian action) are located on the circle at Eu- 
clidian time Ti, each with a 3-momentum conserving 5- 
function; (iii) one integrates over all vertex positions t 
and all propagator 3-momenta p; (iv) a factor of (—1)"^ 
is included where vp is the winding number carried by 
fermions in the diagram, with an additional (—1) for 
each closed fermion loop in the parent Feynman diagram; 
(v) symmetry factors are computed as in Feynman dia- 
grams, with the caveat that two propagators connecting 
the same two vertices do not warrant a symmetry factor 
when their length differs by n/3; (vi) the winding num- 
ber about compact Euclidian time, weighted by particle 
charge, is the order of the graph; all graphs of order p 
are included in a fugacity expansion to order . For ex- 
ample, a dimer loop with v = 1 contributes to order 
since the dimer has particle number 2. 

III. INTERACTIONS AND 62 

To include 2-particle interactions it is convenient to 
represent the interaction not in terms of a potential, but 
by s-channel dimer exchange, a technique introduced in 
[8]. The advantage is that the dimer — with a disper- 
sion relation constructed to exactly reproduce the two 
particle phase shift 5 — has a separable contact interac- 
tion with the fermions. This phase shift is assumed to 
be given, either directly from scattering data, or previ- 
ously calculated from a potential model. That one can 
take this simplifying approach is due to the fact that the 
virial coefficients depend on interactions only through the 
5-matrix [TO]. 

We focus on the case of s-wave scattering, comment- 
ing below on its generalization to other partial waves. 
Consider the Minkowski spacetime Lagrangian 

C = i)\idt + VV2Af)7/' + (l)^K(t) + i^V^o-aV' + h.c. (6) 

where K is a function of the Galilean invariant operator 
D = {idt + V"^ /AM). It is important to recognize that 
this is not a conventional effective field theory (EFT); 
in an EFT, one would perform a low energy expansion 
and express K as a, polynomial in D, as was done to 
subleading order in [8 and to leading order in [U; such 
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FIG. 3: Feynman graphs for dimer mediated two-body scat- 
tering, and the integral equation relevant for three-body scat- 
tering (dashed line — fully dressed dimer, solid line = fermion; 
gray = K'^). 
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FIG. 5: The three propagator subgraph (in black) contribut- 
ing to 63 and corresponding to the expression L appearing in 
eqs. ( [lO|lT| . 



across the cut. Combining the cut and pole contribu- 
tions and performing the P integration we immediately 
recover the well-known result [Til [121 



FIG. 4: Chronograph expansion for the virial coefficients 
bn- Dashed line: fully dressed dimer propagator, shaded box: 
summed three-body interaction from Fig. [S] For 63, bosons 
or multiple fermion species would require introduction of a 
trimer field as well. 



an expansion of K in powers of D corresponds directly to 
the effective range expansion of pcot 5{p). However, here 
we consider X to be a more general function of D (e.g, 
nonlocal) chosen so that the full dimer Green's function 
given by the sum in Fig. |3] results in the exact 2- fermion 
scattering amplitude: 



An 



M kcotS{k) + V^ ' 



(7) 



where k'^ = [M{E + 2fj.) - P^/A + ie], E and P being 
the total energy and momentum of the fermion pair. In 
Fig. [3] the geometric sum of loop diagrams creates the 
correct 2-fermion cut appearing as the V— fc^ term in 
the amplitude; the loops are linearly divergent, and the 
divergence is absorbed into a constant counterterm in 
so that the renormalized operator is = K — const. 
(see, for example [21 [3]). This theory is valid beyond the 
radius of convergence of the effective range expansion, up 
to energies where inelastic processes set in, such as pion 
production in the case where the fermions are nucleons. 

The chronographs in this theory for ^1,2,3 are shown 
in Fig. |4| The second coefficient 62_gets the contribution 
.(1) ^ 



'\/A\f2 computed in eq. (|4j) from free fermions, 
as well as the dimer contribution computed from dQ/d^: 



u(2) 



A3 

4 



(27r)3 



dE 
2^' 



M 



(8) 



where G^^ and k are given in eq. (j?). The contour inte- 
gration picks up contributions botliirom poles and from 
the cut along the positive real E axis from \/ —k'^ in eq. 
([7]). If the theory has bound states with binding energy 
e„, then &^ has poles at E^ = (e„-|-p2/4M-2^) and the 
integrand in brackets has poles a.t E = En with residue 
2. To compute the contribution from the cut, one substi- 
tutes eq. til for G^^ accounting for \/—k^ flipping sign 



u(2) _ 



1 



'dk^m,-^^ 

dk 



(9) 



This analysis can be extended to other partial waves by 
introducing new dimer fields with appropriate couplings 
to fermions. 



IV. COMPUTING 63 

At third order we need to compute the new chrono- 
graphs for n shown in Fig. |4j the first graph yielding the 
free fermion contribution, ~ l/9\/3 as computed in 

eq. ^ . The second graph for 63'^^ sums all dimer- fermion 
and three fermion interactions, barring additional three- 
body forces; for a single species of fermion renormaliz- 
ability does not require three-body forces and we shall 
ignore them here; however, in principle a trimer could be 
introduced to generate fundamental three-body forces. 
This diagram can be expressed in a loop expansion as 



L" 



(10) 



where L corresponds to the subdiagram in Fig. [5] with 

ii2ki,34k2 = Gain - Ti,P/3 - ki) 

X V{t3 - T2, 2P/3 + ki)go(T4 - T3, P/3 + ki + ka) 



(11) 



Go being the fermion chronograph propagator and V be- 
ing the dimer propagator computed from eqs. (5[7), and 
we include a symmetry factor of 1/n, a spin factor of 
Trcrl" ~ 2 from the vertices and a factor of (—1) due to 
i>F = 3; there is no overall fermion sign from the parent 
Feynman diagram. In this expression, a product of L's 
corresponds to the integral 



-^12141,56142 



dr-idTA 



d3q 



^12ki,34qi34q,56k2 (12) 



with an integral over the center of mass momentum P 
implied in the trace; this P integral is gaussian and con- 
tributes a factor of (3\/3/A3). Then the n*'' term in eq. 
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FIG. 6: Results for In/n, including error estimate from nu- 
merical integration for n = 3, . . . , 8. The solid line is the large- 
n fit to the last five data points by the function (co/n+ci/n^), 
as described in the text. 



(10) corresponds to the (n + l)-loop contribution to the 
last diagram in Fig. |4j 

At this point in order to be concrete we will narrow our 

(3) 

focus to computing 63 for unitary fermions, for which 
pcotiJ = identically. From eqs. (5j7) we find the dimer 
propagator X'(T,q) = -6'(r)^167r/Af3 e-'''^^ ''^'^ / ^ . In 
this case the expansion takes the form 



1,(3) ^ 33/2 

3 2-^ ri 



(13) 



where /„ is given by the integral 



(167r 



V2 



n 



n 



(2^)-" / fr4^^^i^(detA)-i 



■J g-k„-A„6kt, 



(14) 



which is positive and independent of (3 and M, A being 
the n X n matrix 



A, 



Sat 



4r2a+2 - T2a+1 + T2a - 4T2a-l 



Sa. 



T2a — T2a-1 



b+1- 



O b 



(15) 



with Sab is the Kronecker (5-function with indices defined 



modulo n, so that 5n+a,b = Sa,n+b = Sab for 1 < a, 6 < n, 
and Jj- represents a 2n-dimensional time ordered integral 
over T,j with < ti < T2 < . . . < r2„ < 1 and T2n+j = 
1 + Tj for j = 1,2. 

The integral /„ can be performed analytically for n = 
1,2 with the result h = 2/{3V5Tr), h = 8/(9V37r); for 
n = 3, . . . , 8 we have computed the integrals numeri- 
cally. We find that /„ is apparently a smooth function 
of n for large n and we perform a large-n extrapola- 
tion, fitting our results for n — 4, . . . , 8 to the func- 
tion /„ ~ (co + ci/n) finding cq = 0.1955 ± 0.0011, 
ci = -0.1943 ± 0.0044 (with highly correlated errors); 
the results are shown in Fig. [6] We estimate our errors 
by also fitting to the function (cq + ci/n -\- C2/n^), and 
varying the range of points used in the fit, finding a very 
stable result. Our final result for the interacting contri- 



bution to 63 is h 



(3) 



-0.3573 ± 0.0005, or for the fuh 



answer, 63 = —0.29315 ± 0.0005, to be compared with 
the recent computation 63 — —0.29095295 by Liu, Hu 
and Drummond [1 , which involved summing over en- 
ergy levels for three unitary fermions in a harmonic trap 
|13) . It is remarkable that an expansion and extrapola- 
tion of chronographs is able to arrive at such a precise 
number for what is essentially a nonperturbative system; 
however, we do not have an explanation for why our pro- 
cedure leads to an estimated error for 63 (±0.0005) which 
is significantly smaller than our discrepancy with Liu et 
al. (0.0021). ' 

To our knowledge, this is the first time 63 was calcu- 
lated by analytical means for a strongly interacting sys- 
tem, and it suggests that the graphical techniques pre- 
sented here for the virial expansion may prove powerful 
for applications to other systems as well. 
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